/*
 * 	This file is a modification of phylip's drawtree.c. It is also part 
 *  of MVD_GUI.
 *
 *  MVD_GUI is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MVD_GUI is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with MVD_GUI.  If not, see <http://www.gnu.org/licenses/>.
 *  This version copyright 2010 Desmond Schmidt
 */
/* version 3.6. (c) Copyright 1993-2004 by the University of Washington.
   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, Andrew Keeffe,
   and Dan Fineman.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <syslog.h>
#include <assert.h>
#include <string.h>
#include "phylip.h"

int spp;
naym *nayme; /* names of species */
/** 
 * Count the number of nodes in a ring, return the total number of 
 * nodes excluding the one passed into the function (siblings)     
 */
int count_sibs (node *p)
{ 
  	node *q;
  	int return_int = 0;
  
  	if (p->tip) 
	{
   		syslog(LOG_ERR,"The function count_sibs called on a tip.  This is a bug.\n");
    	exit (-1);
  	}
  	q = p->next;
  	while (q != p) 
	{
    	if (q == NULL) 
		{
      		syslog(LOG_ERR,"A loop of nodes was not closed.\n");
      		exit (-1);
    	} 
		else 
		{
      		return_int++;
      		q = q->next;
    	}
  	}
  	return return_int;
}
/**
 * Random number generator -- slow but machine independent
 * This is a multiplicative congruential 32-bit generator
 * x(t+1) = 1664525 * x(t) mod 2^32, one that passes the
 * Coveyou-Macpherson and Lehmer tests, see Knuth ACP vol. 2
 * We here implement it representing each integer in base-64
 * notation -- i.e. as an array of 6 six-bit chunks   
 */
double randum(longer seed)
{ 
  	int i, j, k, sum;
  	longer mult, newseed;
  	double x;

  	mult[0] = 13;   /* these four statements set the multiplier */
  	mult[1] = 24;   /* -- they are its "digits" in a base-64    */
  	mult[2] = 22;   /*    notation: 1664525 = 6*64^3+22*64^2    */
  	mult[3] = 6;    /*                         +24*64+13        */
  	for (i = 0; i <= 5; i++)
    	newseed[i] = 0;
	/* do the multiplication piecewise */
  	for (i = 0; i <= 5; i++) 
	{  
    	sum = newseed[i];
    	k = i;
    	if (i > 3)
      		k = 3;
    	for (j = 0; j <= k; j++)
      		sum += mult[j] * seed[i - j];
    	newseed[i] = sum;
    	for (j = i; j <= 4; j++) 
		{
      		newseed[j + 1] += newseed[j] / 64;
      		newseed[j] &= 63;
    	}
  	}
  	memcpy(seed, newseed, sizeof(longer));        /* new seed replaces old one */
  	seed[5] &= 3;          /* from the new seed, get a floating point fraction */
  	x = 0.0;
  	for (i = 0; i <= 5; i++)
    	x = x / 64.0 + seed[i];
  	x /= 4.0;
  	return x;
}
/** 
 * Randomize input order of species.
 * randomly permute array enterorder 
 */
void randumize(longer seed, int *enterorder)
{ 
  	int i, j, k;
  
  	for (i = 0; i < spp; i++) 
	{
		j = (int)(randum(seed)*(i+1));
		k = enterorder[j];
		enterorder[j] = enterorder[i];
		enterorder[i] = k;
  	}
}
/* input random number seed */
void initseed(int *inseed, int *inseed0, longer seed)
{ 
  	int i;

  	assert(inseed);
  	assert(inseed0);
  	*inseed0 = *inseed = 1;
  	for (i = 0; i <= 5; i++)
    	seed[i] = 0;
  	i = 0;
  	do 
	{
    	seed[i] = *inseed & 63;
    	*inseed /= 64;
    	i++;
  	} while (*inseed != 0);
}
/** 
 * Hook together two nodes 
 */
void hookup(node *p, node *q)
{ 
	assert(p != NULL);
	assert(q != NULL);
	p->back = q;
	q->back = p;
} 
/** This loop reads a name from treefile and stores it in *str.
 * Returns the length of the name string. str must be at
 * least MAXNCH bytes, but no effort is made to null-terminate
 * the string. Underscores and newlines are converted to spaces.
 * Characters beyond MAXNCH are discarded. */
int take_name_from_tree(char *ch, char *str, char *treefile, int *i)
{
	int name_length = 0;
	do 
	{
    	if ((*ch) == '_')
      		(*ch) = ' ';
    	if ( name_length < MAXNCH )
      		str[name_length++] = (*ch);
    	if (eoln_str(treefile,*i)) 
      		scan_eoln_str(treefile,i);
    	(*ch) = gettc_str(treefile,i);
    	if (*ch == '\n')
      		*ch = ' ';
  	} 
	while ( strchr(":,)[;", *ch) == NULL );
  	return name_length;
}
/** Recursive procedure adds nodes to user-defined tree
 * This is the main (new) tree-reading procedure 
 */
void addelement(node **p, node *q, char *ch, int *parens, char *treefile,
    int *i, pointarray treenode, boolean *goteof, boolean *first, pointarray nodep,
    int *nextnode, int *ntips, boolean *haslengths, node **grbg,
    initptr initnode, boolean unifok, int maxnodes)
{
  	node *pfirst;
  	int j, len = 0, nodei = 0;
  	boolean notlast;
  	char str[MAXNCH+1];
  	node *r;
  	int furs = 0;

  	if ((*ch) == '(') 
	{
    	(*nextnode)++;          /* get ready to use new interior node */
    	nodei = *nextnode;      /* do what needs to be done at bottom */
    	if ( maxnodes != -1 && nodei > maxnodes) 
		{
			syslog(LOG_ERR,"ERROR in input tree file: Attempting to allocate too\n" 
			"many nodes. This is usually caused by a unifurcation.\n" 
			"To use this tree with this program  use Retree to read\n"
			"and write this tree.\n");
      		exit(-1);
    	}
		/* do what needs to be done at bottom */
    	(*initnode)(p, grbg, q, len, nodei, ntips,
			parens, bottom, treenode, nodep, str, ch, treefile,i);
		pfirst = (*p);
    	notlast = true;
    	while (notlast) 
		{  /* loop through immediate descendants */
      		furs++;
      		(*initnode)(&(*p)->next, grbg, q,
				len, nodei, ntips, parens, nonbottom, treenode,
				nodep, str, ch, treefile, i);
			/* ... doing what is done before each */
			r = (*p)->next;
      		getch_str(ch, parens, treefile, i);      /* look for next character */
			/* handle blank names */
      		if((*ch) == ',' || (*ch) == ':')
			{
        		(*i)--;
				*ch = 0;
      		} 
			else if((*ch)==')')
			{
        		(*i)--;
				(*parens)++;
        		*ch = 0;
			}
      		addelement(&(*p)->next->back, (*p)->next, ch, parens, treefile,
        		i, treenode, goteof, first, nodep, nextnode, ntips,
        		haslengths, grbg, initnode, unifok, maxnodes);  
      		/* this does nothing because the hslength case doesn't occur 
			(*initnode)(&r, grbg, q, len, nodei, ntips,
				parens, hslength, treenode, nodep, str, ch, treefile, i);*/
			 /* do what is done after each about length */
      		pfirst->numdesc++;               /* increment number of descendants */
      		*p = r;                         /* make r point back to p */
      		if ((*ch) == ')') 
			{
        		notlast = false;
        		do 
				{
          			getch_str(ch, parens, treefile,i);
        		} 
				while ((*ch) != ',' && (*ch) != ')' &&
           			(*ch) != '[' && (*ch) != ';' && (*ch) != ':');
			}
		}
    	if ( furs <= 1 && !unifok ) 
		{
      		syslog(LOG_ERR,"ERROR in input tree file: A Unifurcation was detetected.\n"
				"To use this tree with this program use retree to read and"
				" write this tree\n");
      		exit(-1);
		}  
    	(*p)->next = pfirst;
    	(*p) = pfirst;
	} 
	else if ((*ch) != ')') 
	{       /* if it's a species name */
    	for (j = 0; j < MAXNCH+1; j++)   /* fill string with nulls */
      		str[j] = '\0';
    	len = take_name_from_tree (ch, str, treefile, i);  /* get the name */
		if ((*ch) == ')')
      		(*parens)--;         /* decrement count of open parentheses */
    	(*initnode)(p, grbg, q, len, nodei, ntips,
			parens, tip, treenode, nodep, str, ch, treefile, i);
		/* do what needs to be done at a tip */
	} 
	else
    	getch_str(ch, parens, treefile,i);
	if (q != NULL)
    	hookup(q, (*p));                    /* now hook up */
  	(*initnode)(p, grbg, q, len, nodei, ntips, 
		parens, iter, treenode, nodep, str, ch, treefile,i);
	/* do what needs to be done to variable iter */
	if ((*ch) == ':')
    	(*initnode)(p, grbg, q, len, nodei, ntips, 
			parens, length, treenode, nodep, str, ch, treefile, i);
	/* do what needs to be done with length */
	else if ((*ch) != ';' && (*ch) != '[')
    	(*initnode)(p, grbg, q, len, nodei, ntips, 
			parens, hsnolength, treenode, nodep, str, ch, treefile, i);
	/* ... or what needs to be done when no length */
  	if ((*ch) == '[')
    	(*initnode)(p, grbg, q, len, nodei, ntips,
			parens, treewt, treenode, nodep, str, ch, treefile, i);
	/* ... for processing a tree weight */
	else if ((*ch) == ';')     /* ... and at end of tree */
    	(*initnode)(p, grbg, q, len, nodei, ntips,
		parens, unittrwt, treenode, nodep, str, ch, treefile, i);
}
/* Print a message and exit when EOF is reached prematurely. */
void EOF_error()
{ 
  	syslog(LOG_ERR,"Unexpected end-of-file.\n");
  	exit(-1);
}
/** Return the next character in file.
 * If EOF is reached, print an error and die.
 * DOS ('\r\n') and Mac ('\r') newlines are returned as a single '\n'. 
 */
char gettc_str(char* file, int *i) 
{ 
	char ch = file[(*i)++];
	if ( ch == 0 )
    	EOF_error();
	if ( ch == '\r' ) 
	{
    	ch = file[(*i)++];
    	if ( ch != '\n' )
      		(*i)--;
    	ch = '\n';
  	}
  	return ch;
}
/** Return the next character in file.
 * If EOF is reached, print an error and die.
 * DOS ('\r\n') and Mac ('\r') newlines are returned as a single '\n'. 
 */
char gettc(FILE* file) 
{ 
	int ch;
	ch=getc(file);
	if ( ch == EOF )
    	EOF_error();
	if ( ch == '\r' ) 
	{
    	ch = getc(file);
    	if ( ch != '\n' )
      		ungetc(ch, file);
    	ch = '\n';
  	}
  	return ch;
}
/* collect garbage on p -- put it on front of garbage list */
void chuck(node **grbg, node *p)
{
	p->back = NULL;
  	p->next = *grbg;
  	*grbg = p;
}
/**
 * Eat everything up to EOF or newline, including newline
 */
void scan_eoln_str( char *f, int *i ) 
{ 
  	char ch;
	while (!eoff_str(f,*i) && !eoln_str(f,*i)) 
    	gettc_str(f,i);
  	if (!eoff_str(f,*i)) 
    	ch = gettc_str(f,i);
}
/**
 * Eat everything up to EOF or newline, including newline
 */
void scan_eoln( FILE *f ) 
{ 
  	char ch;
  	while (!eoff(f) && !eoln(f)) 
    	gettc(f);
  	if (!eoff(f)) 
    	ch = gettc(f);
}
/** Return true iff strinng is at end */
boolean eoff_str(char *f,int start)
{
    return (f[start]==0);
} 
/** Return true iff next getc() is EOF */
boolean eoff(FILE *f)
{
    int ch;
    if (feof(f)) 
      	return true;
    ch = getc(f);
    if (ch == EOF) 
	{
      	ungetc(ch, f);
      	return true;
	}
    ungetc(ch, f);
    return false;
} 
/* Return true iff next getc() is EOL or EOF */
boolean eoln_str(char *f, int start)
{
    register int ch;
    ch = f[start];
    if (ch == 0)
      	return true;
    return ((ch == '\n') || (ch == '\r'));
}
/* Return true iff next getc() is EOL or EOF */
boolean eoln(FILE *f)
{
    register int ch;
    ch = getc(f);
    if (ch == EOF)
      	return true;
    ungetc(ch, f);
    return ((ch == '\n') || (ch == '\r'));
}
/** 
 * Get next non-blank character from a string
 */
void getch_str(char *c, int *parens, char *treefile, int *i)
{ 
  	do 
	{
    	if (eoln_str(treefile,*i)) 
     		scan_eoln_str(treefile,i);
    	(*c) = gettc_str(treefile,i);
		if ((*c) == '\n' || (*c) == '\t')
      		(*c) = ' ';
  	} 
	while ( *c == ' ' && !eoff_str(treefile,*i) );
  	if ((*c) == '(')
    	(*parens)++;
  	if ((*c) == ')')
    	(*parens)--;
} 
/** 
 * Get next non-blank character
 */
void getch(char *c, int *parens, FILE *treefile)
{ 
  	do 
	{
    	if (eoln(treefile)) 
     		scan_eoln(treefile);
    	(*c) = gettc(treefile);
		if ((*c) == '\n' || (*c) == '\t')
      		(*c) = ' ';
  	} 
	while ( *c == ' ' && !eoff(treefile) );
  	if ((*c) == '(')
    	(*parens)++;
  	if ((*c) == ')')
    	(*parens)--;
} 
/**
 * Read a branch length from a treefile 
 */
void processlength( double *valyew, double *divisor, char *ch, 
	boolean *lengthIsNegative, char *treefile, int *i, int *parens )
{
	int digit, ordzero, exponent, exponentIsNegative;
  	boolean pointread, hasExponent;

  	ordzero = '0';
  	*lengthIsNegative = false;
  	pointread = false;
  	hasExponent = false;
  	exponentIsNegative = -1; // 3 states:  -1 = unassigned, 1 = true, 0 = false
  	exponent = 0;
  	*valyew = 0.0;
  	*divisor = 1.0;
  	getch_str( ch, parens, treefile,i );
  	if ('+' == *ch)
    	getch_str(ch, parens, treefile,i); // ignore leading '+', because "+1.2345" == "1.2345"
  	else if ('-' == *ch)
    {
      	*lengthIsNegative = true;
      	getch_str(ch, parens, treefile,i);
    }
  	digit = (int)(*ch - ordzero);
  	while ( ((digit <= 9) && (digit >= 0)) || '.' == *ch || '-' == *ch
	  	|| '+' == *ch || 'E' == *ch || 'e' == *ch) 
	{
		if ('.' == *ch)
		{
			if (!pointread)
		  		pointread = true;
			else
		  	{
				syslog(LOG_ERR,"Branch length found with more than one \'.\' in it.\n\n");
				exit(-1);
		  	}
		}
		else if ('+' == *ch)
		{
			if (hasExponent && -1 == exponentIsNegative)
		  		exponentIsNegative = 0; // 3 states:  -1 = unassigned, 1 = true, 0 = false
			else
			{
				syslog(LOG_ERR,"Branch length found with \'+\' in an unexpected place.\n\n");
				exit(-1);
			}
		}
		else if ('-' == *ch)
      	{
			if (hasExponent && -1 == exponentIsNegative)
	  			exponentIsNegative = 1; // 3 states:  -1 = unassigned, 1 = true, 0 = false
			else
		  	{
				syslog(LOG_ERR,"Branch length found with \'-\' in an unexpected place.\n\n");
				exit(-1);
		  	}
      	}
    	else if ('E' == *ch || 'e' == *ch)
      	{
			if (!hasExponent)
	  			hasExponent = true;
			else
	  		{
	    		syslog(LOG_ERR,"Branch length found with more than one \'E\' in it.\n\n");
	    		exit(-1);
	  		}
      	}
    	else 
		{
      		if (!hasExponent)
			{
	  			*valyew = *valyew * 10.0 + digit;
	  			if (pointread)
	    			*divisor *= 10.0;
			}
      		else
				exponent = 10*exponent + digit;
    	}
    	getch_str(ch, parens, treefile,i);
    	digit = (int)(*ch - ordzero);
  	}
  	if (hasExponent)
    {
      	if (exponentIsNegative)
			*divisor *= pow(10.,(double)exponent);
      	else
			*divisor /= pow(10.,(double)exponent);
    }
  	if (*lengthIsNegative)
    	*valyew = -(*valyew);
}
/**
 * Initialization of node pointers, variables 
 */
void setupnode( node *p, int i )
{
	p->next = NULL;
	p->back = NULL;
	p->times_in_tree = (double) i * 1.0;
	p->index = i;
	p->tip = false;
}
/**
 * This and the following are do-it-yourself garbage collectors.
 * Make a new node or pull one off the garbage list 
 */
void gnu( node **grbg, node **p )
{ 
  	if (*grbg != NULL) 
	{
    	*p = *grbg;
    	*grbg = (*grbg)->next;
  	} 
	else
    	*p = (node *)Malloc(sizeof(node));
	(*p)->back = NULL;
	(*p)->next = NULL;
	(*p)->tip = false;
	(*p)->times_in_tree = 0.0;
	(*p)->r = 0.0;
	(*p)->theta = 0.0;
	(*p)->x = NULL;
	(*p)->protx = NULL;        /* for the sake of proml     */
} 

/**
 * Read in user-defined tree and set it up. Eats blank lines and 
 * everything up to the first open paren, then calls the recursive 
 * function addelement, which builds the tree and calls back to 
 * initnode. 
 */
void treeread (char *treefile, node **root, pointarray treenode,
	boolean *goteof, boolean *first, pointarray nodep, 
	int *nextnode, boolean *haslengths, node **grbg, initptr initnode,
	boolean unifok, int maxnodes )
{
	char  ch;
  	int parens = 0;
  	int ntips = 0;
	int i = 0;

  	(*goteof) = false;
  	(*nextnode) = spp;

  	/* eat blank lines */
  	while (eoln_str(treefile,i) && !eoff_str(treefile,i))
    	scan_eoln_str(treefile,&i);

  	if (eoff_str(treefile,i)) 
	{
    	*goteof = true;
    	return;
  	} 
  	getch_str(&ch, &parens, treefile,&i);
  	while (ch != '(') 
	{
    	/* Eat everything in the file (i.e. digits, tabs) until you
       	encounter an open-paren */
    	getch_str(&ch, &parens, treefile,&i);
  	}
  	if (haslengths != NULL)
    	*haslengths = true; 
  	addelement(root, NULL, &ch, &parens, treefile, &i,
         treenode, goteof, first, nodep, nextnode, &ntips,
         haslengths, grbg, initnode, unifok, maxnodes);

  	/* Eat blank lines and end of current line*/
  	do 
	{
    	scan_eoln_str(treefile,&i);
  	}
  	while (eoln_str(treefile,i) && !eoff_str(treefile,i));

  	if (first)
    	*first = false;
  	if (parens != 0) {
    	syslog(LOG_ERR,"ERROR in tree file: unmatched parentheses\n");
    	exit(-1);
  	}
} 
void memerror()
{
  	syslog(LOG_ERR,"Error allocating memory\n");
  	exit(-1);
}
/**
 * Error message if attempt to malloc too little or too much memory 
 * @param x the requested amount
 */
void odd_malloc( int x )
{ 
  	syslog(LOG_ERR,"A function asked for an inappropriate amount of memory\n");
  	exit(-1);
}
/**
 * Wrapper for malloc, allowing error message if too little, 
 * or too much
 */
void *mymalloc( int x )
{ 
  	void *new_block;

  	if ((x <= 0) || (x > TOO_MUCH_MEMORY))
    	odd_malloc(x);
  	new_block = (void*)calloc(1, x);
  	if (!new_block) 
	{
   		memerror();
    	return (void*) new_block;
  	} 
	else
    	return (void*) new_block;
}
/**
 * Read the contents of the tree input and allocate as a string
 * @param treeFile input open treefile
 * @param VAR param will return buffer user must free
 */
void readTreeContents( FILE *treeFile, char **treeContents )
{
	int len = 1024;
	int i = 0;
	*treeContents = malloc( len );
	if ( *treeContents == NULL )
	{
		syslog(LOG_ERR,"Failed to allocate treefile buffer\n");
		exit( -1 );
	}
	while ( !feof(treeFile) )
	{
		char c = getc( treeFile );
		if ( !feof(treeFile) )
		{
			if ( i < len )
				(*treeContents)[i++] = c;
			else
			{
				char *temp = malloc( len+1024 );
				if ( temp == NULL )
				{
					syslog(LOG_ERR,"Failed to reallocate treefile buffer\n");
					exit( -1 );			
				}
				else
				{
					int j;
					for ( j=0;j<len;j++ )
						temp[j] = (*treeContents)[j];
					free( *treeContents );
					*treeContents = temp;
					len += 256;
				}
			}
		}
	}
	(*treeContents)[i] = 0;
}
/** 
 * Original routine didn't work with stdin. Rewrote it to use 
 * a string instead.
 * @param treeContents return the treefile contents as a string
 * @param comma return number of commas detected
*/
int countcomma( char *treeContents, int *comma )
{
	int i,len;
	int lparen = 0;
	boolean ignore = false;
	len = strlen( treeContents );
	*comma = 0;
   	for ( i=0;i<len;i++ )
	{
    	if ( !ignore )
		{
			if (treeContents[i] == ',')
		  		(*comma)++;
			if (treeContents[i] == '(')
		     	lparen++;
			if (treeContents[i] == '[') 
		  		ignore = true;
		}
		else if ( treeContents[i] == ']' )
			ignore = false;
	}
  	return lparen + (*comma);
} 
/**
 * Pre-compute space and allocate memory for nodep
 */
void allocate_nodep( pointarray *nodep, char *treeContents, int *precalc_tips )  
{ 
	int numnodes;      /* returns number commas & (    */
	int numcom = 0;        /* returns number commas */
	numnodes = countcomma( treeContents, &numcom ) + 1;
  	*nodep = (pointarray)Malloc( 2*numnodes*sizeof(node*) );
	/* this will be used in placing the tip nodes in the front 
       region of nodep.  Used for species check?  */
  	(*precalc_tips) = numcom + 1;        
}
